#![allow(non_camel_case_types)]
#![allow(non_snake_case)]
extern crate docify;
use super::bsm::Option;
use statrs::distribution::{Continuous, ContinuousCDF, Normal};
use std::f64::consts::E;
use std::f64::consts::PI;
use std::f64::NAN;
const DS: f64 = 0.001;
const DR: f64 = 0.00001;
const DT: f64 = 1.0 / 365.0;
const DV: f64 = 0.00001;

///Barone-Adesi-Whaley定价模型 所有模型来自于 https://zhuanlan.zhihu.com/p/671852660
///https://download.csdn.net/blog/column/12405537/133420200
#[doc = docify::embed!("src/option_price/baw.rs", test_baw)]
/// 当持有成本(Cost of Carry)大于无风险利率时，美式期权与欧式期权的价值相同，
impl Option {
    pub fn baw_call(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let b = self.r - self.divdend;
        let eu = self.bsm_call(vol);
        if b >= self.r {
            return eu;
        }
        let Sk = Kc(self, vol);
        let N = 2.0 * b / vol.powf(2.0);
        let k = 2.0 * self.r / (vol.powf(2.0) * (1.0 - (-1.0 * self.r * self.T).exp()));
        let d1 =
            ((Sk / self.K).log(E) + (b + vol.powf(2.0) / 2.0) * self.T) / (vol * (self.T.sqrt()));
        let Q2 = (-1.0 * (N - 1.0) + ((N - 1.0).powf(2.0) + 4.0 * k)).sqrt() / 2.0;
        let a2 = (Sk / Q2) * (1.0 - ((b - self.r) * self.T).exp() * norm.cdf(d1));
        if self.S < Sk {
            return eu + a2 * (self.S / Sk).powf(Q2);
        }
        self.S - self.K
    }

    pub fn baw_put(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let b = self.r - self.divdend;
        let eu = self.bsm_put(vol);

        let Sk = Kp(self, vol);
        let N = 2.0 * b / vol.powf(2.0);
        let k = 2.0 * self.r / (vol.powf(2.0) * (1.0 - (-1.0 * self.r * self.T).exp()));
        let d1 =
            ((Sk / self.K).log(E) + (b + vol.powf(2.0) / 2.0) * self.T) / (vol * (self.T.sqrt()));
        let Q1 = (-1.0 * (N - 1.0) - ((N - 1.0).powf(2.0) + 4.0 * k).sqrt()) / 2.0;
        let a1 = -1.0 * (Sk / Q1) * (1.0 - ((b - self.r) * self.T).exp() * norm.cdf(-1.0 * d1));
        if self.S > Sk {
            return eu + a1 * (self.S / Sk).powf(Q1);
        }
        self.K - self.S
    }
    pub fn baw_price(&self, vol: f64) -> f64 {
        if self.is_call {
            return self.baw_call(vol);
        } else {
            return self.baw_put(vol);
        }
    }
    pub fn baw_implied_vol(&self, mp: f64) -> f64 {
        let precise = 1e-3;
        let mut l = Vec::new();
        let mut start = 0.5;
        for n in 0..4000 {
            let price = self.baw_price(start);
            let dif = mp - self.baw_price(start);
            if dif.is_nan() {
                continue;
            }
            if dif.abs() < precise {
                return start;
            }
            l.push((dif.abs(), start));
            let mut vega = self.baw_vega(start);
            if vega == 0. {
                vega = 0.01;
            }
            start = start + precise * 2. * dif / vega;
            if start < 0. {
                start = 0.001;
            }
        }
        l.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
        l[0].1
    }
    pub fn baw_delta(&self, vol: f64) -> f64 {
        let mut op1 = self.clone();
        let mut op2 = self.clone();
        op1.S = self.S - DS;
        op2.S = self.S + DS;
        let p1 = op1.baw_price(vol);
        let p2 = op2.baw_price(vol);
        (p2 - p1) / (2.0 * DS)
    }
    pub fn baw_gamma(&self, vol: f64) -> f64 {
        let mut op1 = self.clone();
        let mut op2 = self.clone();
        op1.S = self.S - DS;
        op2.S = self.S + DS;

        let p1 = op1.baw_price(vol);
        let p2 = op2.baw_price(vol);
        let p0 = self.baw_price(vol);
        (p2 - 2.0 * p0 + p1) / (DS.powf(2.0))
    }

    pub fn baw_vega(&self, vol: f64) -> f64 {
        let p1 = self.baw_price(vol - DV);
        let p2 = self.baw_price(vol + DV);
        (p2 - p1) / (2.0 * DV)
    }
    pub fn baw_theta(&self, vol: f64) -> f64 {
        let mut op1 = self.clone();
        op1.T = self.T - DT;
        let p1 = op1.baw_price(vol);
        let p2 = self.baw_price(vol);
        p1 - p2
    }
    pub fn baw_rho(&self, vol: f64) -> f64 {
        let mut op1 = self.clone();
        let mut op2 = self.clone();
        op1.r = self.r - DR;
        op2.r = self.r + DR;
        let p1 = op1.baw_price(vol);
        let p2 = op2.baw_price(vol);
        (p2 - p1) / (2.0 * DR)
    }
}

pub fn Kc(op: &Option, vol: f64) -> f64 {
    let b = op.r - op.divdend;
    let N = 2.0 * b / vol.powf(2.0);
    let m = 2.0 * op.r / vol.powf(2.0);
    let q2u = (-1.0 * (N - 1.0) + ((N - 1.0).powf(2.0) + 4.0 * m).sqrt()) / 2.0;
    let su = op.K / (1.0 - 1.0 / q2u);
    let h2 = -1.0 * (b * op.T + 2.0 * vol * (op.T).sqrt()) * op.K / (su - op.K);
    let mut Si = op.K + (su - op.K) * (1.0 - h2.exp());
    let k = 2.0 * op.r / (vol.powf(2.0) * (1.0 - (-1.0 * op.r * op.T).exp()));
    let mut d1 = ((Si / op.K).log(E) + (b + vol.powf(2.0) / 2.0) * op.T) / (vol * (op.T).sqrt());
    let Q2 = (-1.0 * (N - 1.0) + ((N - 1.0).powf(2.0) + 4.0 * k).sqrt()) / 2.0;
    let mut LHS = Si - op.K;
    let mut op2 = op.clone();
    op2.S = Si;
    let europ = &op2.bsm_call(vol);
    let norm = Normal::standard();
    let mut RHS = europ + (1.0 - ((b - op.r) * op.T).exp() * norm.cdf(d1)) * Si / Q2;
    let mut bi = ((b - op.r) * op.T).exp() * norm.cdf(d1) * (1.0 - 1.0 / Q2)
        + (1.0 - ((b - op.r) * op.T).exp() * norm.pdf(d1) / (vol * (op.T).sqrt())) / Q2;

    while (LHS - RHS).abs() / op.K > 0.001 {
        Si = (op.K + RHS - bi * Si) / (1.0 - bi);
        d1 = ((Si / op.K).log(E) + (b + vol.powf(2.0) / 2.0) * op.T) / (vol * (op.T).sqrt());
        LHS = Si - op.K;
        op2.S = Si;
        RHS = &op2.bsm_call(vol) + (1.0 - ((b - op.r) * op.T).exp() * norm.cdf(d1)) * Si / Q2;
        bi = ((b - op.r) * op.T).exp() * norm.cdf(d1) * (1.0 - 1.0 / Q2)
            + (1.0 - ((b - op.r) * op.T).exp() * norm.cdf(d1) / (vol * (op.T).sqrt())) / Q2;
    }
    Si
}
pub fn Kp(op: &Option, vol: f64) -> f64 {
    let b = op.r - op.divdend;
    let N = 2.0 * b / vol.powf(2.0);
    let m = 2.0 * op.r / vol.powf(2.0);
    let q1u = (-1.0 * (N - 1.0) - ((N - 1.0).powf(2.0) + 4.0 * m).sqrt()) / 2.0;
    let su = op.K / (1.0 - 1.0 / q1u);
    let h1 = (b * op.T - 2.0 * vol * (op.T).sqrt()) * op.K / (op.K - su);
    let mut Si = su + (op.K - su) * h1.exp();
    let k = 2.0 * op.r / (vol.powf(2.0) * (1.0 - (-1.0 * op.r * op.T).exp()));
    let mut d1 = ((Si / op.K).log(E) + (b + vol.powf(2.0) / 2.0) * op.T) / (vol * (op.T).sqrt());
    let Q1 = (-1.0 * (N - 1.0) - ((N - 1.0).powf(2.0) + 4.0 * k).sqrt()) / 2.0;
    let mut LHS = op.K - Si;

    let mut op2 = op.clone();
    op2.S = Si;
    let mut europ = op2.bsm_put(vol);
    let norm = Normal::standard();
    let mut RHS = europ - (1.0 - ((b - op.r) * op.T).exp() * norm.cdf(-1.0 * d1)) * Si / Q1;

    let mut bi = -1.0 * ((b - op.r) * op.T).exp() * norm.cdf(-1.0 * d1) * (1.0 - 1.0 / Q1)
        - (1.0 + ((b - op.r) * op.T).exp() * norm.pdf(-1.0 * d1) / (vol * (op.T).sqrt())) / Q1;

    while (LHS - RHS).abs() / op.K > 0.001 {
        Si = (op.K - RHS + bi * Si) / (1.0 + bi);
        d1 = ((Si / op.K).log(E) + (b + vol.powf(2.0) / 2.0) * op.T) / (vol * (op.T).sqrt());

        LHS = op.K - Si;
        op2.S = Si;
        europ = op2.bsm_put(vol);
        RHS = europ - (1.0 - ((b - op.r) * op.T).exp() * norm.cdf(-1.0 * d1)) * Si / Q1;
        bi = ((b - op.r) * op.T).exp() * norm.cdf(-1.0 * d1) * (1.0 - 1.0 / Q1)
            - (1.0 + ((b - op.r) * op.T).exp() * norm.cdf(-1.0 * d1) / (vol * (op.T).sqrt())) / Q1;
    }
    Si
}

#[test]
#[docify::export]
fn test_baw() {
    // let op = Option {
    //     is_call: true,
    //     S: 220.11, // base asset price
    //     K: 220.0,  // exacuate price
    //     T: 10.0 / 365.0,
    //     divdend: 0.44 / 100.0, //divdent rate
    //     r: 4.07 / 100.0,       //risk free rate
    // };

    // let vol = 28.9 / 100.0;
    // let mp = 4.4;
    // let price = op.baw_price(vol);
    // let theta = op.baw_theta(vol);
    // let vega = op.baw_vega(vol);
    // let rho = op.baw_rho(vol);
    // let gamma = op.baw_gamma(vol);
    // let delta = op.baw_delta(vol);
    // let imply = op.baw_implied_vol(mp);
    // dbg!("baw==", imply, price, delta, gamma, theta, vega, rho);
    let op2 = Option {
        is_call: false,
        S: 73110.,  // base asset price
        K: 75000.0, // exacuate price
        T: 0.12,
        divdend: 0., //divdent rate
        r: 0.013,    //risk free rate
    };
    let mp = 2600.;
    let iv = op2.baw_implied_vol(mp);
    dbg!(iv);
}
